Análisis
Elio
Analizando datos speedy
(luego de la reunión del 5 de agosto de 2016)
Cambiando, decidimos dejar un poco de lado el geopotencial y ver más función corriente, u y T.
Variables: geop, streamf, T, U, V. Cosas para ver: * Estado medio * SD por latitud, mes * Mapas de las ondas * R^2 de las ondas * Cortes longitudinales de cada onda en distintas latitudes
Comparación con NCEP (pasar a la misma resolución usando interpolación bilineal)
Comparación con NCEP:
NCEP no brinda información de función corriente, para comparar uso geopotencial, U y T.
# Leo speedy y hago anomalías mensuales
m_speedy <- fread("SPEEDY/todo_m_speedy.csv")
m_speedy[, month := factor(month, levels = month.abb)]## lon lat lev month gh u v
## 1: 0.00 -87.159 925 Ene 736.8700 -0.05454388 -0.3241713
## 2: 3.75 -87.159 925 Ene 738.2792 -0.02213817 -0.3285298
## 3: 7.50 -87.159 925 Ene 739.8143 0.01149593 -0.3288868
## 4: 11.25 -87.159 925 Ene 741.4750 0.04583307 -0.3251453
## 5: 15.00 -87.159 925 Ene 743.2585 0.08036167 -0.3172835
## ---
## 221180: 341.25 -1.856 30 Dic 23976.2846 -2.46443088 0.3448354
## 221181: 345.00 -1.856 30 Dic 23975.1531 -2.47587266 0.4439990
## 221182: 348.75 -1.856 30 Dic 23972.9128 -2.67345830 0.5403290
## 221183: 352.50 -1.856 30 Dic 23970.0830 -3.03758931 0.4697203
## 221184: 356.25 -1.856 30 Dic 23967.2616 -3.41034904 0.2759294
## psi temp
## 1: -0.1354184 267.7331
## 2: -0.1314101 267.6461
## 3: -0.1278538 267.5543
## 4: -0.1247234 267.4581
## 5: -0.1219879 267.3583
## ---
## 221180: 15.9369396 206.9459
## 221181: 16.0817268 207.1822
## 221182: 16.2614006 207.4649
## 221183: 16.4383895 207.8276
## 221184: 16.5591864 208.0658# setindex(speedy, lon, lat, lev, month, year)
vars <- c("gh", "u", "v", "psi", "temp")
avars <- paste0("a", vars)
allvars <- c(vars, avars)
sp_levs <- m_speedy[, unique(lev)]
# Anomalías mensuales
m_speedy[, c(avars) := lapply(.SD, FUN = anom),
by = .(lat, lev, month),
.SDcols = vars]## lon lat lev month gh u v
## 1: 0.00 -87.159 925 Ene 736.8700 -0.05454388 -0.3241713
## 2: 3.75 -87.159 925 Ene 738.2792 -0.02213817 -0.3285298
## 3: 7.50 -87.159 925 Ene 739.8143 0.01149593 -0.3288868
## 4: 11.25 -87.159 925 Ene 741.4750 0.04583307 -0.3251453
## 5: 15.00 -87.159 925 Ene 743.2585 0.08036167 -0.3172835
## ---
## 221180: 341.25 -1.856 30 Dic 23976.2846 -2.46443088 0.3448354
## 221181: 345.00 -1.856 30 Dic 23975.1531 -2.47587266 0.4439990
## 221182: 348.75 -1.856 30 Dic 23972.9128 -2.67345830 0.5403290
## 221183: 352.50 -1.856 30 Dic 23970.0830 -3.03758931 0.4697203
## 221184: 356.25 -1.856 30 Dic 23967.2616 -3.41034904 0.2759294
## psi temp agh au av
## 1: -0.1354184 267.7331 0.7799904 0.06318588 -0.36023388
## 2: -0.1314101 267.6461 2.1891375 0.09559159 -0.36459240
## 3: -0.1278538 267.5543 3.7242286 0.12922569 -0.36494939
## 4: -0.1247234 267.4581 5.3849729 0.16356284 -0.36120783
## 5: -0.1219879 267.3583 7.1684832 0.19809143 -0.35334602
## ---
## 221180: 15.9369396 206.9459 5.0669393 1.14545083 0.14958942
## 221181: 16.0817268 207.1822 3.9354289 1.13400905 0.24875303
## 221182: 16.2614006 207.4649 1.6951294 0.93642341 0.34508308
## 221183: 16.4383895 207.8276 -1.1346883 0.57229240 0.27447440
## 221184: 16.5591864 208.0658 -3.9561076 0.19953267 0.08068342
## apsi atemp
## 1: -0.0008850379 0.25983430
## 2: 0.0031231956 0.17283174
## 3: 0.0066795316 0.08099112
## 4: 0.0098099039 -0.01514434
## 5: 0.0125453975 -0.11492359
## ---
## 221180: 2.6544129519 -0.83302260
## 221181: 2.7992001999 -0.59670922
## 221182: 2.9788739988 -0.31396132
## 221183: 3.1558628865 0.04871619
## 221184: 3.2766597577 0.28689470setcolorder(m_speedy, c("lon", "lat", "lev", "month", allvars))
setorder(m_speedy, month, -lev, lat, lon)m_ncep <- fread("NCEP/todo_m_ncep.csv")
m_ncep <- m_ncep[lat <= 0]
m_ncep[, month := factor(month, levels = month.abb)]## lon lat lev month gh temp u
## 1: 0.00 -87.159 1000 Ene 3.194028e-01 -3.746012 -1.28944888
## 2: 3.75 -87.159 1000 Ene 2.009247e+00 -3.783916 -0.98486770
## 3: 7.50 -87.159 1000 Ene 3.678018e+00 -3.820472 -0.68162480
## 4: 11.25 -87.159 1000 Ene 5.273846e+00 -3.858975 -0.37856168
## 5: 15.00 -87.159 1000 Ene 6.822608e+00 -3.896709 -0.07922632
## ---
## 470012: 341.25 -1.856 10 Dic 3.092088e+04 -42.738416 -7.69329709
## 470013: 345.00 -1.856 10 Dic 3.092209e+04 -42.620110 -7.63892995
## 470014: 348.75 -1.856 10 Dic 3.092438e+04 -42.696332 -7.64588785
## 470015: 352.50 -1.856 10 Dic 3.092421e+04 -42.711142 -7.58927508
## 470016: 356.25 -1.856 10 Dic 3.092315e+04 -42.988538 -7.17830892
## v
## 1: -3.37433997
## 2: -3.36687496
## 3: -3.33007242
## 4: -3.25966232
## 5: -3.16194629
## ---
## 470012: -0.17117253
## 470013: -0.02299725
## 470014: -0.08035137
## 470015: -0.24902995
## 470016: -0.65007687m_ncep[, psi := NA]## lon lat lev month gh temp u
## 1: 0.00 -87.159 1000 Ene 3.194028e-01 -3.746012 -1.28944888
## 2: 3.75 -87.159 1000 Ene 2.009247e+00 -3.783916 -0.98486770
## 3: 7.50 -87.159 1000 Ene 3.678018e+00 -3.820472 -0.68162480
## 4: 11.25 -87.159 1000 Ene 5.273846e+00 -3.858975 -0.37856168
## 5: 15.00 -87.159 1000 Ene 6.822608e+00 -3.896709 -0.07922632
## ---
## 470012: 341.25 -1.856 10 Dic 3.092088e+04 -42.738416 -7.69329709
## 470013: 345.00 -1.856 10 Dic 3.092209e+04 -42.620110 -7.63892995
## 470014: 348.75 -1.856 10 Dic 3.092438e+04 -42.696332 -7.64588785
## 470015: 352.50 -1.856 10 Dic 3.092421e+04 -42.711142 -7.58927508
## 470016: 356.25 -1.856 10 Dic 3.092315e+04 -42.988538 -7.17830892
## v psi
## 1: -3.37433997 NA
## 2: -3.36687496 NA
## 3: -3.33007242 NA
## 4: -3.25966232 NA
## 5: -3.16194629 NA
## ---
## 470012: -0.17117253 NA
## 470013: -0.02299725 NA
## 470014: -0.08035137 NA
## 470015: -0.24902995 NA
## 470016: -0.65007687 NAm_ncep[, temp := temp + 273.16]## lon lat lev month gh temp u
## 1: 0.00 -87.159 1000 Ene 3.194028e-01 269.4140 -1.28944888
## 2: 3.75 -87.159 1000 Ene 2.009247e+00 269.3761 -0.98486770
## 3: 7.50 -87.159 1000 Ene 3.678018e+00 269.3395 -0.68162480
## 4: 11.25 -87.159 1000 Ene 5.273846e+00 269.3010 -0.37856168
## 5: 15.00 -87.159 1000 Ene 6.822608e+00 269.2633 -0.07922632
## ---
## 470012: 341.25 -1.856 10 Dic 3.092088e+04 230.4216 -7.69329709
## 470013: 345.00 -1.856 10 Dic 3.092209e+04 230.5399 -7.63892995
## 470014: 348.75 -1.856 10 Dic 3.092438e+04 230.4637 -7.64588785
## 470015: 352.50 -1.856 10 Dic 3.092421e+04 230.4489 -7.58927508
## 470016: 356.25 -1.856 10 Dic 3.092315e+04 230.1715 -7.17830892
## v psi
## 1: -3.37433997 NA
## 2: -3.36687496 NA
## 3: -3.33007242 NA
## 4: -3.25966232 NA
## 5: -3.16194629 NA
## ---
## 470012: -0.17117253 NA
## 470013: -0.02299725 NA
## 470014: -0.08035137 NA
## 470015: -0.24902995 NA
## 470016: -0.65007687 NAm_ncep[, c(avars) := lapply(.SD, FUN = anom), by = .(lat, lev, month),
.SDcols = vars]## lon lat lev month gh temp u
## 1: 0.00 -87.159 1000 Ene 3.194028e-01 269.4140 -1.28944888
## 2: 3.75 -87.159 1000 Ene 2.009247e+00 269.3761 -0.98486770
## 3: 7.50 -87.159 1000 Ene 3.678018e+00 269.3395 -0.68162480
## 4: 11.25 -87.159 1000 Ene 5.273846e+00 269.3010 -0.37856168
## 5: 15.00 -87.159 1000 Ene 6.822608e+00 269.2633 -0.07922632
## ---
## 470012: 341.25 -1.856 10 Dic 3.092088e+04 230.4216 -7.69329709
## 470013: 345.00 -1.856 10 Dic 3.092209e+04 230.5399 -7.63892995
## 470014: 348.75 -1.856 10 Dic 3.092438e+04 230.4637 -7.64588785
## 470015: 352.50 -1.856 10 Dic 3.092421e+04 230.4489 -7.58927508
## 470016: 356.25 -1.856 10 Dic 3.092315e+04 230.1715 -7.17830892
## v psi agh au av apsi atemp
## 1: -3.37433997 NA 10.054246 -0.0262046 -4.17046447 NA -0.05855587
## 2: -3.36687496 NA 11.744090 0.2783766 -4.16299946 NA -0.09645984
## 3: -3.33007242 NA 13.412862 0.5816195 -4.12619692 NA -0.13301612
## 4: -3.25966232 NA 15.008689 0.8846826 -4.05578682 NA -0.17151928
## 5: -3.16194629 NA 16.557452 1.1840180 -3.95807079 NA -0.20925289
## ---
## 470012: -0.17117253 NA 3.675745 2.1369449 -0.04305793 NA 0.19750639
## 470013: -0.02299725 NA 4.888170 2.1913120 0.10511735 NA 0.31581262
## 470014: -0.08035137 NA 7.173449 2.1843541 0.04776323 NA 0.23959038
## 470015: -0.24902995 NA 7.003035 2.2409669 -0.12091535 NA 0.22478026
## 470016: -0.65007687 NA 5.947643 2.6519330 -0.52196226 NA -0.05261632np_levs <- m_ncep[, unique(lev)]
setcolorder(m_ncep, c("lon", "lat", "lev", "month", allvars))
setorder(m_ncep, month, -lev, lat, lon)m_speedy.long <- melt(m_speedy, id.vars = c("lon", "lat", "lev", "month"))
m_speedy.interp <- m_speedy.long[,
approx(x = lev, y = value,
xout = np_levs),
by = .(lon, lat, month, variable)]
colnames(m_speedy.interp)[5:6] <- c("lev", "sp")
m_speedy <- m_speedy.interpCorrelación lineal para cada nivel, mes y variable, separando en parte total y parte asimétrica.
m_ncep.long <- melt(m_ncep, id.vars = c("lon", "lat", "lev", "month"), value.name = "nc")
sp_nc <- merge(m_speedy.interp, m_ncep.long)
remove(m_speedy.interp)
cors <- sp_nc[lev %in% sp_levs, .(Correlacion = cor(sp, nc)), by = .(lev, month, variable)]
levs <- unique(cors$lev)
corsa <- cors[grepl("a", variable)]$Correlacion
cors <- cors[!grepl("a", variable)]
cors[, Asimetrica := corsa]## lev month variable Correlacion Asimetrica
## 1: 30 Ene gh -0.8679965 0.6589105
## 2: 100 Ene gh 0.9841290 0.7028101
## 3: 200 Ene gh 0.9856019 0.7279014
## 4: 300 Ene gh 0.9838896 0.6130320
## 5: 500 Ene gh 0.9801555 0.4849896
## ---
## 476: 300 Dic temp 0.9751213 0.5669539
## 477: 500 Dic temp 0.9929736 0.5298177
## 478: 700 Dic temp 0.9819914 0.3370901
## 479: 850 Dic temp 0.9843473 0.6985642
## 480: 925 Dic temp 0.9895365 0.8309316colnames(cors) <- c("lev", "month", "variable", "Total", "Asimetrica")
cors <- melt(cors, id.vars = c("lev", "month", "variable"), variable.name = "Parte",
value.name = "Correlacion")ggplot(cors, aes(lev, Correlacion, color = Parte)) +
geom_hline(yintercept = 0, color = "gray45") +
geom_line() + facet_grid(month~variable) +
scale_color_brewer(type = "qual", name = "Variable") +
coord_flip() + scale_x_continuous(trans = "reverselog", breaks = levs) +
scale_y_continuous(limits = c(-1, 1), minor_breaks = NULL) +
xlab("Nivel") + theme_bw() +
theme(legend.position = "bottom", legend.key.height = unit(5, "points")) +
ggtitle("Correlación lineal entre Speedy y NCEP ")
## Note: no visible binding for global variable 'xend'
## Note: no visible binding for global variable 'yend'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'R'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'
## Note: no visible binding for global variable 'key.row'
## Note: no visible binding for global variable 'key.col'
## Note: no visible binding for global variable 'label.row'
## Note: no visible binding for global variable 'label.col'interp.levs <- exp(seq(log(925), log(30), length.out = 40))
cors.interp <- cors[!is.na(Correlacion),
interp.dt(as.numeric(month), lev, Correlacion, linear = T,
yo = interp.levs),
by = .(Parte, variable)]## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'z'colnames(cors.interp) <- c("Parte", "variable", "month", "lev", "Correlacion")
cor_plot <- ggplot(cors.interp, aes(month, lev)) +
geom_raster(aes(fill = Correlacion)) +
geom_contour(aes(z = Correlacion, color = ..level..),
binwidth = 0.2) +
facet_wrap(variable~Parte, scales = "free", ncol = 2) +
scale_y_continuous(trans = "reverselog", breaks = levs) +
scale_x_continuous(labels = month.abb, breaks = 1:12) +
scale_fill_distiller(type = "div", palette = "RdBu", direction = 1) +
scale_color_gradient(low = "black", high = "black") +
theme_minimal()
# scale_fill_gradient2(low = "#67001f", high = "#053061")
cor_plot <- direct.label(cor_plot, method = "bottom.pieces")Altura geopotencial (gh) Para el caso del campo total, la correlación del campo es buena (>0.8)en casi todos los niveles y meses, excepto en 30hPa durante verano donde los campos ¡está anticorrelacionados! La parte asimétrica zonal muestra valores menores, indicando que gran parte de la correlación del campo total se debe a la capacidad del modelo de reproducir el gradiente meridional. Sin embargo, se siguen obteniendo correlaciones >0.6 en casi todos los niveles y meses. Se observa un mínimo relativo en 500hPa donde en se tienen correlaciones menores durante casi todo el año y uno en niveles altos centrado en septiembre, donde la correlación llega a ser nula.
Viento zonal (U) Las correlaciones con el campo total son >=0.8 en todo el año y todos los niveles, sin embargo, la parte asimétrica muestra correlaciones mucho más baja con un máximo de ~0.6 en 925hPa. Esto indica que el modelo resuelve correctamente la estructura media del Jet, pero no sus variaciones zonales.
Viento meridional (V) Los campos de correlación son prácticamente idénticos entre parte total y parte asimétrica. Ésta muestra un patrón de bajas correlaciones en general, con anticorrelaciones en niveles altos (entre 300 y 200 hPa) durante todo el año
Temperatura (T) La correlación con el campo total muestra una estructura similar que la altura geopotencial, con una excelente correlación en todos los meses para niveles mayores a 200hPa, pero anticorrelacionado en niveles altos en todos los meses salvo en invierno. La parte asimétrica muestra correlaciones bajas en todos los niveles salvo en 925hPa.
Parte simétrica
Altura geopotencial
sp_nc_mean <- sp_nc[, lapply(.SD, mean), by = .(lat, lev, month, variable),
.SDcols = c("sp", "nc")]
sp_nc_mean.long <- melt(sp_nc_mean, id.vars = c("lat", "lev", "month", "variable"), variable.name = "modelo")
ggplot(sp_nc_mean.long[variable == "gh"], aes(lat, value/1000)) +
geom_line(aes(color = modelo, group = interaction(modelo, as.factor(lev)),
linetype = as.factor(lev))) +
scale_linetype_manual(values = rep(c(1,2), length.out = 17)) +
guides(linetype = FALSE) +
geom_text(data = sp_nc_mean.long[month == "Dic" & lat %~% -50 & modelo == "sp" & variable == "gh"], aes(label = lev), nudge_y = 0) +
facet_grid(~month) +
theme_elio +
scale_color_brewer(type = "qual", palette = 6, labels = c("Speedy", "NCEP"),
name = "Modelo") +
scale_x_reverse() +
ylab("Media zonal de altura geopotencial (kmgp)") + xlab("Latitud")
interp.levs <- exp(seq(log(1000), log(10), length.out = 60))
sp_nc_corte <- sp_nc_mean.long[!is.na(value),
interp.dt(lat, lev, value, linear = T,
yo = interp.levs, nx = 60),
by = .(month, variable, modelo)]
colnames(sp_nc_corte)[4:6] <- c("lat", "lev", "value")
sp_nc_corte.wide <- data.table::dcast(sp_nc_corte, month + variable + lat + lev ~ modelo)
corte_latlev_zonal <- function(var, lineas_width, title, scale = "seq") {
M <- sp_nc_corte[variable == var, max(value, na.rm = T)]
m <- sp_nc_corte[variable == var, min(value, na.rm = T)]
m <- floor(m/lineas_width)*lineas_width
M <- ceiling(M/lineas_width)*lineas_width
g <- ggplot(sp_nc_corte.wide[variable == var], aes(lat, lev))
# lineas_width <- 5000
if (scale == "div") {
side <- max(-m, M)
lineas = seq(-side, side, by = lineas_width)
div_pal = colorRampPalette(brewer.pal(name = "RdBu", n = 11))
colors <- div_pal(length(lineas) - 1)
colors <- colors[length(colors):1]
inicio <- 0
g <- g + geom_raster(aes(fill = nc)) +
scale_fill_manual(values = colors, drop = FALSE, name = "NCEP")
} else {
lineas = seq(m, M, by = lineas_width)
inicio <- lineas[1]
g <- g + geom_raster(aes(fill = cut_width(nc, lineas_width, boundary = inicio))) +
scale_fill_viridis(name = "NCEP", option = "plasma",
direction = 1, discrete = T, drop = FALSE,
label = c(lineas, lineas[length(lineas)] + lineas_width))
}
g <- g +
geom_contour(aes(z = sp), alpha = 0.8, color = "black", breaks = lineas) +
geom_dl(aes(z = sp, label = ..level..), method = "top.pieces",
stat = "contour", color = "black", breaks = jumpby(lineas, 2)) +
facet_wrap(~month, ncol = 3) +
scale_y_continuous(trans = "reverselog", breaks = levs) +
scale_x_reverse() +
theme_minimal() +
ggtitle(paste0(title, "\n lineas = ", lineas_width)) +
geom_hline(yintercept = 30, color = "gray45", linetype = 2)
}
g <- corte_latlev_zonal(var = "gh", lineas_width = 5000, title = "Altura geopotencial media (mgp)")## Note: no visible binding for global variable 'variable'
## Note: no visible binding for global variable 'value'
## Note: no visible binding for global variable 'variable'
## Note: no visible binding for global variable 'value'
## Note: no visible binding for global variable 'variable'
## Note: no visible binding for global variable 'lat'
## Note: no visible binding for global variable 'lev'
## Note: no visible binding for global variable 'nc'
## Note: no visible binding for global variable 'nc'
## Note: no visible binding for global variable 'sp'
## Note: no visible binding for global variable 'sp'
## Note: no visible binding for global variable '..level..'g
Comparando la parte zonalmente simétrica de la altura geopotencial se hace evidente la razón de la anticorrelación en niveles altos durante verano. En niveles bajos y medios, el gradiente meridional de geopotencial es positivo (más alturas hacia el ecuador) durane todo el año, con una amplitid que aumenta con la altura. En niveles altos, sin embargo, durante verano NCEP muestra que éste se invierte, mostrando mayores alturas en el altas latitudes. Speedy, por su parte, no logra capturar este comportamiento y sigue repitiendo la estructura básica de nieles inferiores. Esto puede deberse a que Speedy tiene una estratósfera en un nivel más alto que NCEP en los meses de verano.
Más allá de esta salvedad, se ve que todos los niveles existe una buena concordancia entre ambos modelos. La excepción es 50 y 70 hPa, donde Speedy sobreestima la altura geopotencial en todos los meses y latitudes considerablemente.
Temperatura
En la temperatura, se observa que la anticorrelación que se observa en 200 hPa tiene una explicación similar a la de la altura geopotencial. En niveles bajos, el gradiente de temperatura zonal es positivo (mayores temperaturas en el ecuador) mientras que en niveles altos, por encima de 100hPa, este gradiente de invierte. La región intermedia, donde el gradiente es nulo se encuentra cerca de los 200hPa. En los meses de verano para Speedy ésta está lieramente más alta que para NCEP, y ese defasaje es lo que da lugar a la anticorrelación.
Esta diferencia puede indicar que Speedy no captura algún proceso estratosférico. ¿Será que lo que se observa es producto del calentamiento por la producción de Ozono? Hay que investigar cómo Speedy trata el proceso de ozonogénesis, ya que podría ser que no tiene en cuenta el calor extra o lo subestima.
Analizando la parte simétrica de la temperatura, es evidente que Speedy subestima la temperatura significativamente, con diferencias de casi 60K en bajas latitudes en 100hPa. Por otro lado, también sbreestima el ciclo anual en altas latitudes; tanto que en verano invierte el gradiente de temperatura, dando origen a la anticorrelación observada. (En otras palabras, es un desastre)
Observando el perfil en altura.. ¿estratósfera? En -90 speedy parece que ve la estratósfera en en 30hPa mientras que NCEP la pone mucho más abajo, en 300hPa en invierno. En trópico y ecuador no se aprecia la inversión estratosférica.
Speedy sobreestima el gradiente de temperatura. ¿Esto significa más inestabilidad?
plotlevs <- c(850, 500, 300, 200, 100, 30)
ggplot(sp_nc_mean.long[variable == "temp" & lev %in% plotlevs], aes(lat, value)) +
geom_line(aes(color = modelo, group = interaction(modelo, as.factor(lev)))) +
facet_grid(lev~month) +
theme_minimal() +
scale_color_brewer(type = "qual", palette = 6, labels = c("Speedy", "NCEP"),
name = "Modelo") +
scale_y_continuous() +
scale_x_reverse() +
ylab("Media zonal de la temperatura (K)") + xlab("Latitud")
g <- corte_latlev_zonal(var = "temp", lineas_width = 10, title = "Temperatura media (K)")
g
lats <- sp_nc_mean.long[, unique(lat)][c(1, 6, 12, 18, 24)]
ggplot(sp_nc_mean.long[variable == "temp" & lat %in% lats], aes(lev, value - 273.16)) +
geom_line(aes(color = modelo, group = interaction(modelo, as.factor(lat)))) +
facet_grid(lat~month) +
coord_flip() +
theme_minimal() +
scale_color_brewer(type = "qual", palette = 6, labels = c("Speedy", "NCEP"),
name = "Modelo") +
scale_x_continuous(trans = "reverselog", breaks = levs) +
scale_y_continuous(breaks = seq(-90, 30, by = 40)) +
xlab("Nivel (hPa)") + ylab("Temperatura (C)")
Viento zonal
# g <- corte_latlev_zonal(var = "u", lineas_width = 10, title = "Viento zonal (m/s)", scale = "div")
#
# div_pal = colorRampPalette(brewer.pal(name = "RdBu", n = 11))
# g + discrete_scale(values = div_pal(20), name = "NCEP", drop = FALSE)
#
#
plot_var = "u"
M <- ceiling(sp_nc_corte[variable == plot_var, max(value, na.rm = T)])
m <- floor(sp_nc_corte[variable == plot_var, min(value, na.rm = T)])
lineas_width <- 10
m <- floor(m/lineas_width)*lineas_width
lineas = seq(m, M, by = lineas_width)
ggplot(sp_nc_corte.wide[variable == plot_var], aes(lat, lev)) +
geom_raster(aes(fill = nc)) +
geom_contour(aes(z = sp, linetype = as.factor(-sign(..level..))),
alpha = 0.8, color = "black", breaks = lineas) +
geom_dl(aes(z = sp, label = ..level..), method = "top.pieces",
stat = "contour", color = "black", breaks = jumpby(lineas, 2)) +
geom_contour(aes(z = nc),
alpha = 0.8, color = "grey33", breaks = lineas) +
geom_dl(aes(z = nc, label = ..level..), method = "bottom.pieces",
stat = "contour", color = "grey33", breaks = jumpby(lineas, 2)) +
scale_linetype_manual(values = c(1, 2, 3)) +
facet_wrap(~month, ncol = 3) +
scale_y_continuous(trans = "reverselog", breaks = levs) +
scale_x_reverse() +
scale_fill_gradient2(name = "NCEP", high = muted("red"), low = muted("blue")) +
theme_minimal() + ggtitle(paste0("Viento zonal (m/s) \n lineas = ", lineas_width))
Por su parte, el viento zonal ambos modelos muestran una estructura similar con el jet subtropical bien definido aunque Speedy lo muestra más intenso y ligeramente corrido hacia el ecuador y en niveles más altos que NCEP durante los meses de verano. Durante el invierno, NCEP muestra el jet subpolar, pero Speedy no logra representarlo por falta de niveles verticales, aunque la estructura por debajo del jet se muestra similar aunque subestimando la magnitud.
Speedy tampoco muestra los vientos del este en niveles altos que dominan bajas latitudes en invierno y casi todas las latitudes en verano. Tiene un mínimo y valores negativos, pero muy subestimados.
ggplot(sp_nc_mean.long[variable == "u" & lev %in% plotlevs], aes(lat, value)) +
geom_line(aes(color = modelo, group = interaction(modelo, as.factor(lev)))) +
facet_grid(lev~month) +
theme_minimal() +
scale_color_brewer(type = "qual", palette = 6, labels = c("Speedy", "NCEP"),
name = "Modelo") +
scale_y_continuous(breaks = seq(-10, 60, by = 20)) +
scale_x_reverse() +
ylab("Media zonal del viento zonal (m/s)") + xlab("Latitud") +
theme_elio
Viento meridional
plot_var = "v"
M <- ceiling(sp_nc_corte[variable == plot_var, max(value, na.rm = T)])
m <- floor(sp_nc_corte[variable == plot_var, min(value, na.rm = T)])
lineas_width <- .5
m <- floor(m/lineas_width)*lineas_width
lineas = seq(lineas_width, M, by = lineas_width)
ggplot(sp_nc_corte.wide[variable == plot_var], aes(lat, lev)) +
geom_raster(aes(fill = nc)) +
geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
alpha = 0.9, color = "black", breaks = fill_neg(lineas)) +
geom_dl(aes(z = sp, label = ..level..), method = "top.pieces",
stat = "contour", color = "black", breaks = jumpby(fill_neg(lineas), 2)) +
geom_contour(aes(z = nc),
alpha = 0.6, color = "grey33", breaks = fill_neg(lineas)) +
geom_dl(aes(z = nc, label = ..level..), method = "bottom.pieces",
stat = "contour", color = "grey33", breaks = jumpby(fill_neg(lineas), 2)) +
scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
facet_wrap(~month, ncol = 3) +
scale_y_continuous(trans = "reverselog", breaks = levs) +
scale_x_reverse() +
scale_fill_gradient2(name = "NCEP", high = muted("red"), low = muted("blue")) +
theme_elio + ggtitle(paste0("Viento zonal (m/s) \n lineas = ", lineas_width))
plot_var = "v"
M <- ceiling(sp_nc_corte[variable == plot_var, max(value, na.rm = T)])
m <- floor(sp_nc_corte[variable == plot_var, min(value, na.rm = T)])
lineas_width <- .5
m <- floor(m/lineas_width)*lineas_width
lineas = seq(lineas_width, M, by = lineas_width)
ggplot(sp_nc_corte.wide[variable == plot_var & lev >= 500], aes(lat, lev)) +
geom_raster(aes(fill = nc)) +
geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
alpha = 0.9, color = "black", breaks = fill_neg(lineas)) +
geom_dl(aes(z = sp, label = ..level..), method = "top.pieces",
stat = "contour", color = "black", breaks = jumpby(fill_neg(lineas), 2)) +
geom_contour(aes(z = nc),
alpha = 0.6, color = "grey33", breaks = fill_neg(lineas)) +
geom_dl(aes(z = nc, label = ..level..), method = "bottom.pieces",
stat = "contour", color = "grey33", breaks = jumpby(fill_neg(lineas), 2)) +
scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
facet_wrap(~month, ncol = 3) +
scale_y_continuous(trans = "reverselog", breaks = levs) +
scale_x_reverse() +
scale_fill_gradient2(name = "NCEP", high = muted("red"), low = muted("blue")) +
theme_elio + ggtitle(paste0("Viento zonal (m/s) \n lineas = ", lineas_width))
En el viento meridional, Speedy parece capturar la estructura general de la circulación pero subestima las magnitudes. En bajas latitudes, hay un dipolo entre niveles bajos y altos que alterna entre invierno y verano; se trata de la parte convergente en superficie y divergente en altura de la ITCZ, que se mueve hacia el hemisferio de verano.
En altas latitudes, en superficie hay máximos de viento del sur debido a los vientos catabáticos de la antártida. Nuevamente, Speedy captura esta característica, pero con una considerable subestimación.
ggplot(sp_nc_mean.long[variable == "v" & lev %in% plotlevs], aes(lat, value)) +
geom_line(aes(color = modelo, group = interaction(modelo, as.factor(lev)))) +
facet_grid(lev~month, scales = "free_y") +
theme_elio +
scale_color_brewer(type = "qual", palette = 6, labels = c("Speedy", "NCEP"),
name = "Modelo") +
scale_y_continuous() +
scale_x_reverse() +
ylab("Media zonal de la temperatura (K)") + xlab("Latitud")
Parte Antisimétrica
Altura geopotencial
Campos por estación y nivel:
sp_nc[, estacion := asign_season(month)]## lon lat month variable lev sp nc
## 1: 0.00 -87.159 Ene gh 10 NA 31955.3946220
## 2: 0.00 -87.159 Ene gh 20 NA 27002.9561155
## 3: 0.00 -87.159 Ene gh 30 23577.7435547 24162.6416817
## 4: 0.00 -87.159 Ene gh 50 21357.7593099 20622.8639138
## 5: 0.00 -87.159 Ene gh 70 19137.7750651 18308.5150735
## ---
## 4700156: 356.25 -1.856 Dic atemp 600 -0.4605941 -0.5462326
## 4700157: 356.25 -1.856 Dic atemp 700 -1.5416182 -0.5649992
## 4700158: 356.25 -1.856 Dic atemp 850 -1.3644017 -0.4677048
## 4700159: 356.25 -1.856 Dic atemp 925 -1.8917754 -0.4087908
## 4700160: 356.25 -1.856 Dic atemp 1000 NA -0.8071430
## estacion
## 1: Verano
## 2: Verano
## 3: Verano
## 4: Verano
## 5: Verano
## ---
## 4700156: Verano
## 4700157: Verano
## 4700158: Verano
## 4700159: Verano
## 4700160: Veranog_data <- rep_lon(sp_nc[variable == "agh", .(sp = mean(sp), nc = mean(nc)),
by = .(lat, lon, lev, estacion)])## Note: no visible binding for global variable 'lon'
## Note: no visible binding for global variable 'lon'maplevs <- plotlevs[plotlevs != 850]
m_data <- g_data[lev %in% maplevs]
M <- m_data[, (max(c(abs(min(nc)), max(nc))))]
M <- ceiling(M)
# cuts <- seq(-M, M, length.out = 10)
# lineas <- floor(abs(cuts[1] - cuts[2]))
# cuts <- cuts[cuts != 0]
s <- 0.5
linea <- 60
lineas <- seq(linea, M, by = linea)
# m_data <- m_data[lev == 30 & estacion == "Primavera"]
ggplot(m_data, aes(lon, lat)) +
# geom_tile(aes(fill = nc)) +
facet_grid(lev~estacion) +
HS3 +
geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "red4", size = s, alpha = 0.8) +
# geom_dl(aes(z = sp, label = ..level..), method = "top.pieces", stat = "contour", breaks = fill_neg(jumpby(lineas, 2)), color = "red4") +
geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "blue4", size = s, alpha = 0.8) +
# geom_dl(aes(z = nc, label = ..level..), method = "bottom.pieces", stat = "contour", breaks = fill_neg(jumpby(lineas, 2)), color = "blue4") +
scale_fill_gradient2(high = muted("red"), low = muted("blue"),
name = "NCEP") +
scale_linetype_manual(values = c(2, 1, 1)) +
coord_map("stereographic", orientation = c(-90,0, 60),
ylim = c(-90, -20)) +
scale_x_continuous(breaks = seq(0, 330, by = 30)) +
# annotate(geom = "text", x = seq(0, 330, by = 30), y = -5,
# label = seq(0, 330, by = 30), color = "gray45") +
# scale_y_continuous(breaks = seq(0+15, -90+15, by = -15)) +
# annotate(geom = "text", x = 0, y = seq(0+15, -90+15, by = -15),
# label = seq(0+15, -90+15, by = -15), color = "gray45") +
theme_elio +
theme(legend.position = "bottom", legend.key.height = unit(5, "points"),
axis.title = element_blank(), axis.text = element_blank()) +
ggtitle(paste0("Z* para NCEP (azul) y SPEEDY (rojo) \n lineas = ", linea))
Corte en -65º
ggplot(g_data[lat %~% -65], aes(lon, lev)) +
geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "red4", size = s, alpha = 0.8) +
geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "blue4", size = s, alpha = 0.8) +
scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
facet_grid(~estacion) +
scale_y_continuous(trans = "reverselog", breaks = levs) +
theme_bw() +
ggtitle(paste0("Z* para NCEP (azul) y SPEEDY (rojo) en -65º \n lineas = ", linea)) 
Un pequeño experimento no del todo exitoso para visualizar el corte:
ggplot(g_data[lat %~% -65], aes(lon, lev)) +
geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "red4", size = s, alpha = 0.8) +
geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "blue4", size = s, alpha = 0.8) +
scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
facet_grid(~estacion) +
scale_y_continuous(trans = "reverselog", breaks = levs, limits = c(15000, 10)) +
scale_x_continuous(breaks = seq(60, 360, by = 60)) +
theme_bw() +
ggtitle(paste0("Z* para NCEP (azul) y SPEEDY (rojo) en -65º \n lineas = ", linea)) +
coord_polar(start = -300*pi/180, direction = 1) + # annotate(geom = "line", x = c(0, 360), y = rep(925, 2)) +
annotate(geom = "rect", xmax = 360, xmin = 0, ymax = 15000, ymin = 925, fill = "white")
## Note: no visible global function definition for 'scale_transform'Desvío estándar de anomalía de geopotencial.
sp_nc_sd <- sp_nc[, .(sp = sd(sp), nc = sd(nc)),
by = .(lat, month, lev, variable)]
g_data <- sp_nc_sd[lev %in% plotlevs & variable == "agh"]
M <- g_data[, (max(c(abs(min(nc)), max(nc))))]
M <- ceiling(M)
s <- 0.7
linea <- 40
lineas <- seq(linea, M, by = linea)
ggplot(g_data, aes(as.numeric(month), lat)) +
facet_grid(lev~.) +
geom_contour(aes(z = sp), breaks = fill_neg(lineas),
color = "red4", size = s, alpha = 0.8) +
geom_dl(aes(label = ..level.., z = sp), color = "red4",
stat = "contour", method = "bottom.pieces",
breaks = jumpby(fill_neg(lineas), 2)) +
geom_contour(aes(z = nc), breaks = fill_neg(lineas),
color = "blue4", size = s, alpha = 0.8) +
geom_dl(aes(label = ..level.., z = nc), color = "blue4",
stat = "contour", method = "top.pieces",
breaks = jumpby(fill_neg(lineas), 2)) +
theme_elio +
scale_x_continuous(breaks = 1:12, labels = month.abb)
Viento zonal.
g_data <- rep_lon(sp_nc[variable == "au", .(sp = mean(sp), nc = mean(nc)),
by = .(lat, lon, lev, estacion)])
maplevs <- plotlevs[plotlevs != 850]
m_data <- g_data[lev %in% maplevs]
M <- m_data[, (max(c(abs(min(nc)), max(nc))))]
M <- ceiling(M)
s <- 0.3
linea <- 5
lineas <- seq(linea, M, by = linea)
ggplot(m_data, aes(lon, lat)) +
# geom_tile(aes(fill = nc)) +
facet_grid(lev~estacion) +
HS3+
geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "red4", size = s, alpha = 0.8) +
# geom_dl(aes(z = sp, label = ..level..), method = "top.pieces", stat = "contour", breaks = fill_neg(jumpby(lineas, 2)), color = "red4") +
geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "blue4", size = s, alpha = 0.8) +
# geom_dl(aes(z = nc, label = ..level..), method = "bottom.pieces", stat = "contour", breaks = fill_neg(jumpby(lineas, 2)), color = "blue4") +
scale_fill_gradient2(high = muted("red"), low = muted("blue"),
name = "NCEP") +
scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
coord_map("stereographic", orientation = c(-90,0, 60),
ylim = c(-90, 0)) +
scale_x_continuous(breaks = seq(0, 330, by = 30)) +
# annotate(geom = "text", x = seq(0, 330, by = 30), y = -5,
# label = seq(0, 330, by = 30), color = "gray45") +
# scale_y_continuous(breaks = seq(0+15, -90+15, by = -15)) +
# annotate(geom = "text", x = 0, y = seq(0+15, -90+15, by = -15),
# label = seq(0+15, -90+15, by = -15), color = "gray45") +
theme_bw() +
theme(legend.position = "bottom", legend.key.height = unit(5, "points"),
axis.title = element_blank(), axis.text = element_blank()) +
ggtitle(paste0("U* para NCEP (azul) y SPEEDY (rojo) \n lineas = ", linea))
linea <- 2
lineas <- seq(linea, M, by = linea)
ggplot(g_data[lat %~% -65], aes(lon, lev)) +
geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "red4", size = s) +
geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "blue4", size = s) +
scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
facet_grid(~estacion) +
scale_y_continuous(trans = "reverselog", breaks = levs) +
theme_bw() +
ggtitle(paste0("U* para NCEP (azul) y SPEEDY (rojo) en -65º \n lineas = ", linea))
g_data <- sp_nc_sd[lev %in% plotlevs & variable == "au"]
M <- g_data[, (max(c(abs(min(nc)), max(nc))))]
M <- ceiling(M)
s <- 0.7
linea <- 4.5
lineas <- seq(linea, M, by = linea)
ggplot(g_data, aes(as.numeric(month), lat)) +
facet_grid(lev~.) +
geom_contour(aes(z = sp), breaks = fill_neg(lineas),
color = "red4", size = s, alpha = 0.8) +
geom_dl(aes(label = ..level.., z = sp), color = "red4",
stat = "contour", method = "bottom.pieces",
breaks = jumpby(fill_neg(lineas), 2)) +
geom_contour(aes(z = nc), breaks = fill_neg(lineas),
color = "blue4", size = s, alpha = 0.8) +
geom_dl(aes(label = ..level.., z = nc), color = "blue4",
stat = "contour", method = "top.pieces",
breaks = jumpby(fill_neg(lineas), 2)) +
theme_elio +
scale_x_continuous(breaks = 1:12, labels = month.abb)
Temperatura
g_data <- rep_lon(sp_nc[variable == "atemp", .(sp = mean(sp), nc = mean(nc)),
by = .(lat, lon, lev, estacion)])
maplevs <- plotlevs[plotlevs != 850]
m_data <- g_data[lev %in% maplevs]
M <- m_data[, (max(c(abs(min(nc)), max(nc))))]
M <- ceiling(M)
s <- 0.3
linea <- 2
lineas <- seq(linea, M, by = linea)
# m_data <- m_data[lev == 200]
ggplot(m_data, aes(lon, lat)) +
# geom_tile(aes(fill = nc)) +
facet_grid(lev~estacion) +
HS3+
geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "red4", size = s, alpha = 0.8) +
# geom_dl(aes(z = sp, label = ..level..), method = "top.pieces", stat = "contour", breaks = fill_neg(jumpby(lineas, 2)), color = "red4") +
geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "blue4", size = s, alpha = 0.8) +
# geom_dl(aes(z = nc, label = ..level..), method = "bottom.pieces", stat = "contour", breaks = fill_neg(jumpby(lineas, 2)), color = "blue4") +
scale_fill_gradient2(high = muted("red"), low = muted("blue"),
name = "NCEP") +
scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
coord_map("stereographic", orientation = c(-90,0, 60),
ylim = c(-90, 0)) +
scale_x_continuous(breaks = seq(0, 330, by = 30)) +
# annotate(geom = "text", x = seq(0, 330, by = 30), y = -5,
# label = seq(0, 330, by = 30), color = "gray45") +
# scale_y_continuous(breaks = seq(0+15, -90+15, by = -15)) +
# annotate(geom = "text", x = 0, y = seq(0+15, -90+15, by = -15),
# label = seq(0+15, -90+15, by = -15), color = "gray45") +
theme_bw() +
theme(legend.position = "bottom", legend.key.height = unit(5, "points"),
axis.title = element_blank(), axis.text = element_blank()) +
ggtitle(paste0("T* para NCEP (azul) y SPEEDY (rojo) \n lineas = ", linea))
linea <- 2
lineas <- seq(linea, M, by = linea)
ggplot(g_data[lat %~% -65], aes(lon, lev)) +
geom_contour(aes(z = sp, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "red4", size = s, alpha = 0.8) +
geom_contour(aes(z = nc, linetype = as.factor(sign(..level..))),
breaks = fill_neg(lineas), color = "blue4", size = s, alpha = 0.8) +
scale_linetype_manual(values = c(2, 1, 1), guide = "none") +
facet_grid(~estacion) +
scale_y_continuous(trans = "reverselog", breaks = levs) +
theme_bw() +
ggtitle(paste0("T* para NCEP (azul) y SPEEDY (rojo) en -65º \n lineas = ", linea))
g_data <- sp_nc_sd[lev %in% plotlevs & variable == "atemp"]
M <- g_data[, (max(c(abs(min(nc)), max(nc))))]
M <- ceiling(M)
s <- 0.7
linea <- 1
lineas <- seq(linea, M, by = linea)
ggplot(g_data, aes(as.numeric(month), lat)) +
facet_grid(lev~.) +
geom_contour(aes(z = sp), breaks = fill_neg(lineas),
color = "red4", size = s, alpha = 0.8) +
geom_dl(aes(label = ..level.., z = sp), color = "red4",
stat = "contour", method = "bottom.pieces",
breaks = jumpby(fill_neg(lineas), 2)) +
geom_contour(aes(z = nc), breaks = fill_neg(lineas),
color = "blue4", size = s, alpha = 0.8) +
geom_dl(aes(label = ..level.., z = nc), color = "blue4",
stat = "contour", method = "top.pieces",
breaks = jumpby(fill_neg(lineas), 2)) +
theme_elio +
scale_x_continuous(breaks = 1:12, labels = month.abb)
Ondas cuasiestacionarias
colnames(m_ncep.long)[6] <- "value"
colnames(m_speedy.long)[6] <- "value"
m_ncep.long[, modelo := "nc"]## lon lat lev month variable value modelo
## 1: 0.00 -87.159 1000 Ene gh 0.31940276 nc
## 2: 3.75 -87.159 1000 Ene gh 2.00924682 nc
## 3: 7.50 -87.159 1000 Ene gh 3.67801808 nc
## 4: 11.25 -87.159 1000 Ene gh 5.27384569 nc
## 5: 15.00 -87.159 1000 Ene gh 6.82260795 nc
## ---
## 4700156: 341.25 -1.856 10 Dic atemp 0.19750639 nc
## 4700157: 345.00 -1.856 10 Dic atemp 0.31581262 nc
## 4700158: 348.75 -1.856 10 Dic atemp 0.23959038 nc
## 4700159: 352.50 -1.856 10 Dic atemp 0.22478026 nc
## 4700160: 356.25 -1.856 10 Dic atemp -0.05261632 ncm_speedy.long[, modelo := "sp"]## lon lat lev month variable value modelo
## 1: 0.00 -87.159 925 Ene gh 736.87003174 sp
## 2: 3.75 -87.159 925 Ene gh 738.27917887 sp
## 3: 7.50 -87.159 925 Ene gh 739.81427002 sp
## 4: 11.25 -87.159 925 Ene gh 741.47501424 sp
## 5: 15.00 -87.159 925 Ene gh 743.25852458 sp
## ---
## 2211836: 341.25 -1.856 30 Dic atemp -0.83302260 sp
## 2211837: 345.00 -1.856 30 Dic atemp -0.59670922 sp
## 2211838: 348.75 -1.856 30 Dic atemp -0.31396132 sp
## 2211839: 352.50 -1.856 30 Dic atemp 0.04871619 sp
## 2211840: 356.25 -1.856 30 Dic atemp 0.28689470 spsp_nc.long <- rbind(m_ncep.long, m_speedy.long)
remove(m_ncep.long, m_speedy.long)
temp <- sp_nc.long[stringi::stri_sub(variable, 1, 1) == "a"]
qs <- temp[, qs_fit(value, n = 1:4), by = .(lat, lev, month, modelo, variable)]
qs_labels <- c("1" = "QS 1",
"2" = "QS 2",
"3" = "QS 3",
"4" = "QS 4")
mod_labels <- c("nc" = "NCEP",
"sp" = "SPEEDY")R^2
qs_gh <- qs[variable == "agh"]
qs_gh_interp <- qs_gh[!is.na(R.sqr),
interp.dt(lat, lev, R.sqr, linear = T,
yo = interp.levs),
by = .(month, modelo, variable, QS)]
colnames(qs_gh_interp)[5:7] <- c("lat", "lev", "R.sqr")
temp <- qs_gh[!is.na(Amplitud),
interp.dt(lat, lev, Amplitud, linear = T,
yo = interp.levs),
by = .(month, modelo, variable, QS)]
colnames(temp)[5:7] <- c("lat", "lev", "Amplitud")
qs_gh_interp$Amplitud <- temp$Amplitud
g_data <- qs_gh_interp
# g_data <- g_data[lat >= -50 & lat < -35 & lev > 100 & QS == 3]
g_nc <- g_data[modelo == "nc"]
g_sp <- g_data[modelo == "sp"]
lineas_width <- 0.25
lineas <- seq(0, 1, by = lineas_width)
s <- 0.7
ggplot(g_nc, aes(lat, lev)) + scale_y_continuous(trans = "reverselog") +
facet_grid(month ~ QS, labeller = labeller(QS = qs_labels)) +
geom_tile(aes(fill = cut_width(R.sqr, lineas_width, boundary = 0))) +
scale_fill_brewer(palette = "Blues", name = "NCEP") +
geom_contour(aes(z = R.sqr), breaks = lineas,
color = "black", data = g_sp, size = s, alpha = 0.7) +
geom_dl(aes(z = R.sqr, label = ..level..), data = g_sp, color = "black",
stat = "contour", method = list("top.pieces", cex = 0.8),
breaks = jumpby(lineas, 2)) +
theme_elio + ggtitle(paste0("R^2 de cada número de onda por latitud y mes \nNCEP (sombreado) y SPEEDY (contornos) - lineas = ", lineas_width))
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'width'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'width'
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'height'
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'height'Amplitud
s <- 0.7
g_data[lev >= 30 & lev <= 925, MaxAmpl := max(Amplitud), by = .(QS, modelo)]## month modelo variable QS lat lev R.sqr Amplitud
## 1: Ene nc agh 1 -87.159000 1000 0.9773782 23.03031
## 2: Ene nc agh 1 -84.971744 1000 0.9283385 32.59569
## 3: Ene nc agh 1 -82.784487 1000 0.8739121 40.36194
## 4: Ene nc agh 1 -80.597231 1000 0.8079078 44.26125
## 5: Ene nc agh 1 -78.409974 1000 0.7822487 45.99572
## ---
## 230396: Dic sp agh 4 -10.605026 10 NA NA
## 230397: Dic sp agh 4 -8.417769 10 NA NA
## 230398: Dic sp agh 4 -6.230513 10 NA NA
## 230399: Dic sp agh 4 -4.043256 10 NA NA
## 230400: Dic sp agh 4 -1.856000 10 NA NA
## MaxAmpl
## 1: NA
## 2: NA
## 3: NA
## 4: NA
## 5: NA
## ---
## 230396: NA
## 230397: NA
## 230398: NA
## 230399: NA
## 230400: NAg_data[, MaxAmpl := max(MaxAmpl, na.rm = T), by = .(QS, modelo)]## month modelo variable QS lat lev R.sqr Amplitud
## 1: Ene nc agh 1 -87.159000 1000 0.9773782 23.03031
## 2: Ene nc agh 1 -84.971744 1000 0.9283385 32.59569
## 3: Ene nc agh 1 -82.784487 1000 0.8739121 40.36194
## 4: Ene nc agh 1 -80.597231 1000 0.8079078 44.26125
## 5: Ene nc agh 1 -78.409974 1000 0.7822487 45.99572
## ---
## 230396: Dic sp agh 4 -10.605026 10 NA NA
## 230397: Dic sp agh 4 -8.417769 10 NA NA
## 230398: Dic sp agh 4 -6.230513 10 NA NA
## 230399: Dic sp agh 4 -4.043256 10 NA NA
## 230400: Dic sp agh 4 -1.856000 10 NA NA
## MaxAmpl
## 1: 548.81118
## 2: 548.81118
## 3: 548.81118
## 4: 548.81118
## 5: 548.81118
## ---
## 230396: 23.22319
## 230397: 23.22319
## 230398: 23.22319
## 230399: 23.22319
## 230400: 23.22319g_data[, Norm_Ampl := Amplitud/MaxAmpl]## month modelo variable QS lat lev R.sqr Amplitud
## 1: Ene nc agh 1 -87.159000 1000 0.9773782 23.03031
## 2: Ene nc agh 1 -84.971744 1000 0.9283385 32.59569
## 3: Ene nc agh 1 -82.784487 1000 0.8739121 40.36194
## 4: Ene nc agh 1 -80.597231 1000 0.8079078 44.26125
## 5: Ene nc agh 1 -78.409974 1000 0.7822487 45.99572
## ---
## 230396: Dic sp agh 4 -10.605026 10 NA NA
## 230397: Dic sp agh 4 -8.417769 10 NA NA
## 230398: Dic sp agh 4 -6.230513 10 NA NA
## 230399: Dic sp agh 4 -4.043256 10 NA NA
## 230400: Dic sp agh 4 -1.856000 10 NA NA
## MaxAmpl Norm_Ampl
## 1: 548.81118 0.04196399
## 2: 548.81118 0.05939327
## 3: 548.81118 0.07354432
## 4: 548.81118 0.08064932
## 5: 548.81118 0.08380974
## ---
## 230396: 23.22319 NA
## 230397: 23.22319 NA
## 230398: 23.22319 NA
## 230399: 23.22319 NA
## 230400: 23.22319 NAM <- ceiling(g_data[, max(Norm_Ampl, na.rm = T)])
lineas_width <- 0.25
lineas <- c(seq(0, 1, by = lineas_width), 2)
region_3 <- c(latmin = -65, latmax = -40, levmin = 30, levmax = 925)
ggplot(g_data[modelo == "sp"], aes(lat, lev)) +
scale_y_continuous(trans = "reverselog") +
scale_x_reverse() +
facet_grid(month ~ QS, labeller = labeller(QS = qs_labels)) +
geom_raster(data = g_data[modelo == "nc"],
aes(fill = cut(Norm_Ampl, breaks = lineas))) +
geom_contour(aes(z = Norm_Ampl), binwidth = lineas_width,
color = "black", size = s, alpha = 0.7) +
# geom_dl(aes(z = Norm_Ampl, label = ..level..),
# binwidth = lineas_width, stat = "contour", method = "top.pieces") +
annotate(geom = "rect",
xmin = region_3[1], xmax = region_3[2],
ymin = region_3[3], ymax = region_3[4],
fill = NA, color = "black") +
scale_fill_brewer(palette = "Blues", name = "NCEP") +
theme_elio + ggtitle(paste0("Amplitud normalizada para cada número de onda por latitud y mes \nNCEP (sombreado) y SPEEDY (contornos) - lineas = ", lineas_width))
QS3
Nos vamos a quedar con la onda 3 y vamos a analizar la zona entre 65º y 40ºS y entre 925 y 30 hPa
g_data <- qs_gh_interp[QS == 3]
region_3 <- c(latmin = -65, latmax = -40, levmin = 30, levmax = 925)
M <- ceiling(g_data[, max(Amplitud, na.rm = T)])
lineas_width <- 10
lineas <- seq(0, M, by = lineas_width)
ggplot(g_data[modelo == "sp"], aes(lat, lev)) +
scale_y_continuous(trans = "reverselog") +
scale_x_reverse() +
facet_grid(month ~ .) +
geom_raster(data = g_data[modelo == "nc"],
aes(fill = cut(Amplitud, breaks = lineas))) +
geom_contour(aes(z = Amplitud), binwidth = lineas_width,
color = "black", size = s, alpha = 0.7) +
# geom_dl(aes(z = Amplitud, label = ..level..),
# binwidth = lineas_width, stat = "contour", method = "top.pieces") +
annotate(geom = "rect",
xmin = region_3[1], xmax = region_3[2],
ymin = region_3[3], ymax = region_3[4],
fill = NA, color = "black") +
scale_fill_brewer(palette = "Blues", name = "NCEP") +
theme_elio + ggtitle(paste0("Amplitud de QS3 por latitud y mes \nNCEP (sombreado) y SPEEDY (contornos) - lineas = ", lineas_width))
En esa caja puedo tomar el valor promedio (que sería equivalente a la integral, ya que el área se mantiene igual) o el máximo. El ciclo anual de ambas variables tanto para NCEP como para SPEEDY se muestra en la siguiente figura:
qs3_gh <- qs_gh[QS == 3, ]
qs3_gh[, region := (lat >= region_3[1] & lat <= region_3[2] &
lev >= region_3[3] & lev <= region_3[4])]## lat lev month modelo variable Amplitud Fase R.sqr
## 1: -87.159 1000 Ene nc agh 1.122137 -1.8936592 0.002320359
## 2: -83.479 1000 Ene nc agh 8.521985 -1.6079967 0.042457977
## 3: -79.777 1000 Ene nc agh 20.398952 -1.2987663 0.155878052
## 4: -76.070 1000 Ene nc agh 23.695890 -1.0897797 0.203066383
## 5: -72.362 1000 Ene nc agh 11.180591 -1.0035023 0.080182527
## ---
## 7196: -16.700 30 Dic sp agh 4.062847 0.7509012 0.067805701
## 7197: -12.989 30 Dic sp agh 5.123547 1.2430367 0.130525363
## 7198: -9.278 30 Dic sp agh 6.235376 1.3842273 0.261522205
## 7199: -5.567 30 Dic sp agh 7.580177 1.3940119 0.363977906
## 7200: -1.856 30 Dic sp agh 9.331789 1.4057353 0.362829617
## QS region
## 1: 3 FALSE
## 2: 3 FALSE
## 3: 3 FALSE
## 4: 3 FALSE
## 5: 3 FALSE
## ---
## 7196: 3 FALSE
## 7197: 3 FALSE
## 7198: 3 FALSE
## 7199: 3 FALSE
## 7200: 3 FALSEg_data <- qs3_gh[region == T,
.(Mean_Ampl = mean(Amplitud), Max_Ampl = max(Amplitud)),
by = .(month, modelo)]
g_data[, Area_equiv := Mean_Ampl/Max_Ampl]## month modelo Mean_Ampl Max_Ampl Area_equiv
## 1: Ene nc 20.255716 40.10602 0.5050542
## 2: Feb nc 27.311712 54.67475 0.4995306
## 3: Mar nc 22.386954 43.29004 0.5171387
## 4: Abr nc 21.981106 40.32700 0.5450717
## 5: May nc 17.497137 29.57796 0.5915600
## 6: Jun nc 21.332002 36.20745 0.5891606
## 7: Jul nc 25.663798 40.75579 0.6296970
## 8: Ago nc 27.150506 41.26745 0.6579157
## 9: Sep nc 22.962262 36.94286 0.6215616
## 10: Oct nc 15.597001 27.31040 0.5711011
## 11: Nov nc 6.285683 13.24662 0.4745122
## 12: Dic nc 8.606986 18.03878 0.4771379
## 13: Ene sp 12.438160 24.64187 0.5047572
## 14: Feb sp 13.158030 24.46971 0.5377273
## 15: Mar sp 10.394263 18.50236 0.5617804
## 16: Abr sp 12.507457 21.28683 0.5875677
## 17: May sp 17.952138 30.12335 0.5959543
## 18: Jun sp 25.581265 41.36971 0.6183574
## 19: Jul sp 17.528378 27.78134 0.6309408
## 20: Ago sp 15.358337 20.92206 0.7340740
## 21: Sep sp 18.168745 31.32297 0.5800454
## 22: Oct sp 17.962713 31.59355 0.5685563
## 23: Nov sp 22.239114 38.35022 0.5798953
## 24: Dic sp 17.636349 31.80561 0.5545044
## month modelo Mean_Ampl Max_Ampl Area_equivg <- ggplot(g_data, aes(as.numeric(month), Area_equiv)) +
geom_line(aes(color = modelo)) +
scale_color_brewer(type = "qual", palette = 6, labels = c("NCEP", "SPEEDY"),
name = "Modelo", direction = -1) +
scale_x_continuous(breaks = 1:12, labels = month.abb, name = "Mes") +
scale_y_continuous(limits = c(NA, 1)) +
ylab("Amplitud media / Amplitud máxima") +
theme_elio
g_data <- melt(g_data[, !"Area_equiv", with = F],
id.vars = c("modelo", "month"),
variable.name = "Tipo", value.name = "Amplitud")
ggplot(g_data, aes(as.numeric(month), Amplitud)) +
geom_line(aes(color = modelo, linetype = Tipo)) +
scale_x_continuous(breaks = 1:12, labels = month.abb, name = "Mes") +
scale_color_brewer(type = "qual", palette = 6, labels = mod_labels,
name = "Modelo", direction = -1) +
theme_elio
Se observa que no hay mucha correlación entre los modelos. Mientras que NCEP tiene un ciclo semianual con máximos en febrero y agosto, Speedy tiene los máximos en junio y noviembre. Casi una anticorrelación perfecta.
Otra medida de posible relevancia es la relación entre la amplitud máxima y la media. Ésto da una idea de cuan concentrado está la anomalía de geopotencial. Valores cercanos a 1 implican que la misma está distribuida de forma pareja en todo el área de estudio, mientras que valores menores implican una distribución más concentrada.
g
Esta variable muestra un comportamiento muy parecido en ambos modelos, con un máximo en el invierno y un mínimo en verano, indicando que el patrón de QS3 se encuentra mas localizado en verano que en invierno.
Una cosa que hay que tener cuidado es cómo caracterizar la amplitud de la onda 3 en toda la región mediante un número para hacer una serie temporal. La elección más naive es la media, pero ésto sólo es útil si la variable tiene una distribución simétrica y unimodal. Una variable como la amplitud está acotado inferiormente por el cero, por lo que la primera propiedad no se cumple a priori (aunque puede ser aproximada si la amplitud es grande).
La funciones de distribución para cada mes estimadas a partir de los datos muestran que en varios meses tampoco se cumple la segunda propiedad.
ggplot(qs3_gh[region == T, ], aes(Amplitud)) +
geom_density(aes(color = modelo)) + facet_grid(month~.) +
scale_color_brewer(type = "qual", palette = 6, labels = mod_labels,
name = "Modelo", direction = -1) +
xlab("Densidad") +
theme_elio
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'ymax'
## Note: no visible binding for global variable 'ymin'En particular, agosto y noviembre muestran distribuciones multimodales en speedy y ncep respectivamente. Esto quizás podría poner en duda la validez de usar la media como representativo de la actividad de la onda 3 en toda la región.
Serie temporal onda 3
Independientemente de esto, armo una serie temporal mensual con la amplitud de la QS3 para el período 1985-2015 usando datos de NCEP.
# load("NCEP/todo_ncep.rda")
#
# ncep_r3 <- ncep[lat >= region_3[1] & lat <= region_3[2] &
# lev >= region_3[3] & lev <= region_3[4]]
load("NCEP/ncep_r3.rda")
ncep_r3[, agh := anom(gh), by = .(lat, lev, month)]## lon lat lev time gh temp u
## 1: 0.00 -64.942 925 1985-01-01 531.9597 -4.524887 -1.828449
## 2: 3.75 -64.942 925 1985-01-01 521.6587 -4.454431 -2.684299
## 3: 7.50 -64.942 925 1985-01-01 512.3531 -4.373037 -3.581535
## 4: 11.25 -64.942 925 1985-01-01 505.1609 -4.290367 -4.418997
## 5: 15.00 -64.942 925 1985-01-01 500.0706 -4.237465 -5.100862
## ---
## 3499772: 341.25 -42.678 30 2015-12-01 23937.4658 -51.288152 3.380353
## 3499773: 345.00 -42.678 30 2015-12-01 23935.6581 -51.204924 3.560365
## 3499774: 348.75 -42.678 30 2015-12-01 23935.0648 -51.132407 3.708036
## 3499775: 352.50 -42.678 30 2015-12-01 23934.8835 -50.978878 3.867348
## 3499776: 356.25 -42.678 30 2015-12-01 23934.4127 -50.824664 4.091359
## v year month agh
## 1: 4.5704375 1985 1 16.855892
## 2: 4.0168368 1985 1 6.554825
## 3: 3.1929850 1985 1 -2.750730
## 4: 2.2034509 1985 1 -9.942983
## 5: 1.2306730 1985 1 -15.033255
## ---
## 3499772: 0.8035898 2015 12 -4.778157
## 3499773: 0.6716072 2015 12 -6.585878
## 3499774: 0.4757195 2015 12 -7.179177
## 3499775: 0.2354143 2015 12 -7.360489
## 3499776: 0.1725060 2015 12 -7.831237qs3 <- ncep_r3[, qs_fit(agh, n = 3), by = .(lat, lev, time)]
qs3[, month := as.numeric(stringi::stri_sub(time, 6, 7))]## lat lev time Amplitud Fase R.sqr QS month
## 1: -64.942 925 1985-01-01 11.803611 2.074606069 0.115521622 3 1
## 2: -61.232 925 1985-01-01 17.083049 2.266734636 0.114262180 3 1
## 3: -57.521 925 1985-01-01 22.758725 2.442488409 0.142396274 3 1
## 4: -53.810 925 1985-01-01 24.283161 2.554906552 0.155590425 3 1
## 5: -50.099 925 1985-01-01 21.845354 2.589264691 0.154252904 3 1
## ---
## 36452: -57.521 30 2015-12-01 25.537137 -0.337952733 0.012160818 3 12
## 36453: -53.810 30 2015-12-01 20.496370 -0.387828235 0.011492796 3 12
## 36454: -50.099 30 2015-12-01 15.693286 -0.406047472 0.011419940 3 12
## 36455: -46.389 30 2015-12-01 10.640019 -0.271296775 0.010503649 3 12
## 36456: -42.678 30 2015-12-01 5.997953 0.006903003 0.008447028 3 12qs3[, month := factor(month, levels = c(12, 1:11), ordered = T)]## lat lev time Amplitud Fase R.sqr QS month
## 1: -64.942 925 1985-01-01 11.803611 2.074606069 0.115521622 3 1
## 2: -61.232 925 1985-01-01 17.083049 2.266734636 0.114262180 3 1
## 3: -57.521 925 1985-01-01 22.758725 2.442488409 0.142396274 3 1
## 4: -53.810 925 1985-01-01 24.283161 2.554906552 0.155590425 3 1
## 5: -50.099 925 1985-01-01 21.845354 2.589264691 0.154252904 3 1
## ---
## 36452: -57.521 30 2015-12-01 25.537137 -0.337952733 0.012160818 3 12
## 36453: -53.810 30 2015-12-01 20.496370 -0.387828235 0.011492796 3 12
## 36454: -50.099 30 2015-12-01 15.693286 -0.406047472 0.011419940 3 12
## 36455: -46.389 30 2015-12-01 10.640019 -0.271296775 0.010503649 3 12
## 36456: -42.678 30 2015-12-01 5.997953 0.006903003 0.008447028 3 12qs3[, year := as.numeric(stringi::stri_sub(time, 1, 4))]## lat lev time Amplitud Fase R.sqr QS month
## 1: -64.942 925 1985-01-01 11.803611 2.074606069 0.115521622 3 1
## 2: -61.232 925 1985-01-01 17.083049 2.266734636 0.114262180 3 1
## 3: -57.521 925 1985-01-01 22.758725 2.442488409 0.142396274 3 1
## 4: -53.810 925 1985-01-01 24.283161 2.554906552 0.155590425 3 1
## 5: -50.099 925 1985-01-01 21.845354 2.589264691 0.154252904 3 1
## ---
## 36452: -57.521 30 2015-12-01 25.537137 -0.337952733 0.012160818 3 12
## 36453: -53.810 30 2015-12-01 20.496370 -0.387828235 0.011492796 3 12
## 36454: -50.099 30 2015-12-01 15.693286 -0.406047472 0.011419940 3 12
## 36455: -46.389 30 2015-12-01 10.640019 -0.271296775 0.010503649 3 12
## 36456: -42.678 30 2015-12-01 5.997953 0.006903003 0.008447028 3 12
## year
## 1: 1985
## 2: 1985
## 3: 1985
## 4: 1985
## 5: 1985
## ---
## 36452: 2015
## 36453: 2015
## 36454: 2015
## 36455: 2015
## 36456: 2015Meses de verano donde activo e inactivo se definen si la amplitud es mayor, menor o que la media (del mes) en 1 desvío estandar (del mes):
qs3_mean <- qs3[, .(Mean_Ampl = mean(Amplitud),
Max_Ampl = max(Amplitud)),
by = .(month, year, time)]
qs3_mean[, `:=`(MeanMean_Ampl = mean(Mean_Ampl),
SDMean_Ampl = sd(Mean_Ampl),
MeanMax_Ampl = mean(Max_Ampl),
SDMax_Ampl = sd(Max_Ampl)),
by = .(month)]## month year time Mean_Ampl Max_Ampl MeanMean_Ampl SDMean_Ampl
## 1: 1 1985 1985-01-01 24.52334 58.86156 28.67268 11.79553
## 2: 2 1985 1985-02-01 30.08895 67.91443 34.54357 12.04654
## 3: 3 1985 1985-03-01 36.82800 64.30854 33.86206 12.32747
## 4: 4 1985 1985-04-01 38.25564 72.74778 37.35963 17.38566
## 5: 5 1985 1985-05-01 21.27393 32.00394 36.03709 16.72060
## ---
## 368: 8 2015 2015-08-01 24.59589 39.53484 50.04477 23.54863
## 369: 9 2015 2015-09-01 69.22008 100.90178 44.63204 24.43700
## 370: 10 2015 2015-10-01 44.67534 64.20312 34.86898 16.19144
## 371: 11 2015 2015-11-01 15.24799 35.95592 34.99821 15.94231
## 372: 12 2015 2015-12-01 15.31068 28.22762 26.72547 15.62461
## MeanMax_Ampl SDMax_Ampl
## 1: 59.79659 22.14225
## 2: 69.93723 23.27483
## 3: 64.56587 23.92974
## 4: 68.57183 28.63779
## 5: 61.30037 23.99550
## ---
## 368: 75.21750 33.08701
## 369: 71.29075 34.84362
## 370: 59.18704 20.95578
## 371: 60.07507 23.82482
## 372: 51.54154 26.40450shift_season <- function(month, year) {
ifelse(month == 12, year + 1, year)
}
qs3_mean[, Año_estacion := shift_season(month, year)]## month year time Mean_Ampl Max_Ampl MeanMean_Ampl SDMean_Ampl
## 1: 1 1985 1985-01-01 24.52334 58.86156 28.67268 11.79553
## 2: 2 1985 1985-02-01 30.08895 67.91443 34.54357 12.04654
## 3: 3 1985 1985-03-01 36.82800 64.30854 33.86206 12.32747
## 4: 4 1985 1985-04-01 38.25564 72.74778 37.35963 17.38566
## 5: 5 1985 1985-05-01 21.27393 32.00394 36.03709 16.72060
## ---
## 368: 8 2015 2015-08-01 24.59589 39.53484 50.04477 23.54863
## 369: 9 2015 2015-09-01 69.22008 100.90178 44.63204 24.43700
## 370: 10 2015 2015-10-01 44.67534 64.20312 34.86898 16.19144
## 371: 11 2015 2015-11-01 15.24799 35.95592 34.99821 15.94231
## 372: 12 2015 2015-12-01 15.31068 28.22762 26.72547 15.62461
## MeanMax_Ampl SDMax_Ampl Año_estacion
## 1: 59.79659 22.14225 1985
## 2: 69.93723 23.27483 1985
## 3: 64.56587 23.92974 1985
## 4: 68.57183 28.63779 1985
## 5: 61.30037 23.99550 1985
## ---
## 368: 75.21750 33.08701 2015
## 369: 71.29075 34.84362 2015
## 370: 59.18704 20.95578 2015
## 371: 60.07507 23.82482 2015
## 372: 51.54154 26.40450 2016qs3_mean[, Anom := ifelse(Mean_Ampl > MeanMean_Ampl + SDMean_Ampl,
"Activo",
ifelse(Mean_Ampl < MeanMean_Ampl - SDMean_Ampl,
"Inactivo", "Normal"))]## month year time Mean_Ampl Max_Ampl MeanMean_Ampl SDMean_Ampl
## 1: 1 1985 1985-01-01 24.52334 58.86156 28.67268 11.79553
## 2: 2 1985 1985-02-01 30.08895 67.91443 34.54357 12.04654
## 3: 3 1985 1985-03-01 36.82800 64.30854 33.86206 12.32747
## 4: 4 1985 1985-04-01 38.25564 72.74778 37.35963 17.38566
## 5: 5 1985 1985-05-01 21.27393 32.00394 36.03709 16.72060
## ---
## 368: 8 2015 2015-08-01 24.59589 39.53484 50.04477 23.54863
## 369: 9 2015 2015-09-01 69.22008 100.90178 44.63204 24.43700
## 370: 10 2015 2015-10-01 44.67534 64.20312 34.86898 16.19144
## 371: 11 2015 2015-11-01 15.24799 35.95592 34.99821 15.94231
## 372: 12 2015 2015-12-01 15.31068 28.22762 26.72547 15.62461
## MeanMax_Ampl SDMax_Ampl Año_estacion Anom
## 1: 59.79659 22.14225 1985 Normal
## 2: 69.93723 23.27483 1985 Normal
## 3: 64.56587 23.92974 1985 Normal
## 4: 68.57183 28.63779 1985 Normal
## 5: 61.30037 23.99550 1985 Normal
## ---
## 368: 75.21750 33.08701 2015 Inactivo
## 369: 71.29075 34.84362 2015 Activo
## 370: 59.18704 20.95578 2015 Normal
## 371: 60.07507 23.82482 2015 Inactivo
## 372: 51.54154 26.40450 2016 Normal# Alternativa: usando terciles.
# qs3_mean[, Anom := cut(Mean_Ampl,
# breaks = quantile(Mean_Ampl, seq(0, 1, length.out = 4)),
# include.lowest = T, labels = c("Inactivo", "Normal", "Activo")),
# by = month]
years <- seq(1984, 2016, by = 1)
lab <- stri_sub(years, 3, 4)
labs <- paste0(lab, "/", shift(lab, type = "lead"))
n <- length(years)
years <- years[-c(1)]
labs <- labs[-c(n)]
# years <- unique(qs3_mean$Año_estacion)
g <- ggplot(qs3_mean,
aes(Año_estacion, Mean_Ampl)) +
geom_col(aes(fill = Anom)) +
geom_line(aes(y = MeanMean_Ampl), linetype = 2, color = "gray4") +
scale_x_continuous(breaks = jumpby(years, 2),
labels = jumpby(labs, 2),
name = "Año de fin de estación") +
scale_fill_brewer(palette = "Set1") +
ylab("Amplitud media QS3") +
theme_elio +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))
g + facet_grid(month~., labeller = labeller(month = month.abb))
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'width'
## Note: no visible binding for global variable 'x'
## Note: no visible binding for global variable 'width'La gran mayoría de los meses de veranos son normales con 5-6 veranos activo o inactivos.
Se observa que no hay muchos meses con amplitud mucho mayor que la media y que no están agrupados en un mismo verano. Es decir, no hay mucha correlación entre la amplitud de un mes y el siguiente.
Para ver esto último con más exactitud, se puede estimar la función de autocorrelación.
autocor <- acf.sig(qs3_mean$Mean_Ampl, method = "large.lag")
ggplot(autocor, aes(lag, acf)) +
geom_line() + geom_point() +
geom_line(aes(y = sig.cut), linetype = 2) +
geom_line(aes(y = -sig.cut), linetype = 2) +
geom_hline(yintercept = 0) +
scale_x_continuous(limits = c(0, 12*3),
name = "Lag (meses)",
breaks = seq(0, 12*3, by = 6)) +
ylab("Autocorrelación mensual QS3") +
theme_elio
Se ve que la autocorrelación a lag 1 no es significativa (utilizando método large lag), pero además se ve un evidente ciclo anual.
ggplot(qs3_mean, aes(factor(month, labels = month.abb), Mean_Ampl)) +
geom_boxplot() +
xlab("Mes") + ylab("Amplitud media QS3") +
theme_elio
## Note: no visible binding for global variable 'weight'
## Note: no visible binding for global variable 'xmin'
## Note: no visible binding for global variable 'xmax'
## Note: no visible binding for global variable 'y'
## Note: no visible binding for global variable 'size'¿Tiene otros ciclos? Le quito el ciclo anual restando la media anual y hago fourier.
fourier <- as.dt(convert.fft(fft(qs3_mean[, .(Mean_Ampl - MeanMean_Ampl)]$V1)))
ggplot(fourier, aes(freq, ampl)) +
geom_line() +
theme_elio + xlab("Frecuencia") + ylab("Amplitud")
No, no hay ningún ciclo obvio.
qs3 <- merge(qs3, qs3_mean[, .(month, year, Anom)])
#
# g_data <- qs3[, .(Amplitud = mean(Amplitud)), by = .(lat, lev, Anom, month)]
#
# g_data <- g_data[, interp.dt(lat, lev, Amplitud, linear = T,
# yo = interp.levs),
# by = .(Anom, month)]
# colnames(g_data) <- c("Anom", "month", "lat", "lev", "Amplitud")
# # g_data <- g_data[month %in% 1:2]
# ggplot(g_data, aes(lat, lev)) +
# scale_y_continuous(trans = "reverselog") +
# geom_contour(aes(z = Amplitud), binwidth = 10,
# color = "black") +
# geom_dl(aes(z = Amplitud, label = ..level..),
# stat = "contour", binwidth = 20, color = "black",
# method = "top.pieces") +
# facet_grid(month~Anom, labeller = labeller(month = month.abb)) +
# theme_elio
# Esto es la distribución espacial de la amplitud en la región de estudio para cada mes y clasificación. No sé si sirve para mucho. Se ve que en febrero en algunos meses hay más separación que en otros.
ggplot(qs3, aes(Amplitud)) +
geom_density(aes(color = Anom)) +
facet_wrap(~month, ncol = 4, labeller = labeller(month = month.abb)) +
theme_elio
Composiciones de variables.
Ahora la idea es ver los campos de variables que podrían ser forantes y comparar años activos vs. inactivos.
SST
load("NCEP/sst_todo.rda")
## Agrego máscara de océanos
library(maptools)
library(maps)
make_mask <- function(lat, lon) {
seamask <- map("world2", fill=TRUE, col = "transparent", plot = F)
IDs <- sapply(strsplit(seamask$names, ":"), function(x) x[1])
seamask <- map2SpatialPolygons(seamask, IDs = IDs,
proj4string = CRS("+proj=longlat +datum=WGS84"))
points <- SpatialPoints(expand.grid(lon, lat),
proj4string = CRS(proj4string(seamask)))
sea <- is.na(over(points, seamask))
return(sea)
}
lat <- unique(sst$lat)
lon <- unique(sst$lon)
sea <- make_mask(lat, lon)
invisible(sst[, sea := sea])
sst <- merge(sst, qs3_mean[, .(time, Anom, Mean_Ampl)], all = T)
invisible({
sst[, month := as.numeric(stringi::stri_sub(time, 6, 7))]
sst[, month := factor(month, levels = c(12, 1:11), ordered = T)]
sst[sea == F, sst := NA]
sst[, msst := mean(sst, na.rm = T), by = .(lat, month)]
sst[, asst := sst - msst]
})sst_comp <- sst[!is.na(Anom),
.(asst = mean(asst, na.rm = T)),
by = .(lat, lon, Anom, month)]
tmp <- sst_comp[, .(Anom = "Diferencia",
asst = asst[Anom == "Inactivo"] - asst[Anom == "Activo"]),
by = .(lat, lon, month)]
sst_comp <- rbind(sst_comp, tmp)
remove(tmp)
gdata <- sst_comp[Anom %in% c("Activo", "Inactivo")]
g <- ggplot(gdata, aes(lon, lat)) +
geom_contour(aes(z = asst, color = ..level..), binwidth = 2.5) +
world3 +
coord_quickmap() +
scale_color_gradient2(high = muted("red"), low = muted("blue")) +
theme_elioDEF: (líneas = 2.5)
g + facet_grid_paginate(month~Anom, labeller = labeller(month = month.abb),
nrow = 3, ncol = 2,
page = 1)
MAM
g + facet_grid_paginate(month~Anom, labeller = labeller(month = month.abb),
nrow = 3, ncol = 2,
page = 2)
JJA
g + facet_grid_paginate(month~Anom, labeller = labeller(month = month.abb),
nrow = 3, ncol = 2,
page = 3)
SON:
g + facet_grid_paginate(month~Anom, labeller = labeller(month = month.abb),
nrow = 3, ncol = 2,
page = 4)
No se puede ver a ojo mucha diferencia. Restando ambos campos (íneas = 0.5):
lineas <- seq(-3, 3, by = 0.5)
gdata <- sst_comp[Anom == "Diferencia"]
ggplot(gdata, aes(lon, lat)) +
geom_contour(aes(z = asst, color = ..level..), breaks = lineas) +
world3 +
coord_quickmap() +
facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) +
scale_color_gradient2(high = muted("red"), low = muted("blue")) +
theme_elio
En algunos meses se ve que la zona del pacífico ecuatorial tiene grandes diferencias, con nomviembre, diciemre y enero caracterizado por diferencias positivas (años activos con mayor temperatura que inactivos) aunque lo contrario en septiembre y febrero.
Al hacer la clasificación binaria entre activos e inactivos se pierde información. Una alternativa es hacer la regresión de la SST en la amplitud de la onda QS3. Es decir, modelar \(QS3 = m\times SST + b + \epsilon\) para cada punto del globo.
sst[!is.na(Mean_Ampl) & !is.na(asst), sd := sd(asst), by = .(lat, lon, month)]## time lon lat sst sea Anom Mean_Ampl month
## 1: 1984-01-01 0.00 -87.159 NA FALSE NA NA 1
## 2: 1984-01-01 3.75 -87.159 NA FALSE NA NA 1
## 3: 1984-01-01 7.50 -87.159 NA FALSE NA NA 1
## 4: 1984-01-01 11.25 -87.159 NA FALSE NA NA 1
## 5: 1984-01-01 15.00 -87.159 NA FALSE NA NA 1
## ---
## 1769468: 2015-12-01 341.25 87.159 -1.742222 TRUE Normal 15.31068 12
## 1769469: 2015-12-01 345.00 87.159 -1.776820 TRUE Normal 15.31068 12
## 1769470: 2015-12-01 348.75 87.159 -1.790000 TRUE Normal 15.31068 12
## 1769471: 2015-12-01 352.50 87.159 -1.790000 TRUE Normal 15.31068 12
## 1769472: 2015-12-01 356.25 87.159 -1.790000 TRUE Normal 15.31068 12
## msst asst sd
## 1: NaN NA NA
## 2: NaN NA NA
## 3: NaN NA NA
## 4: NaN NA NA
## 5: NaN NA NA
## ---
## 1769468: -1.786385 0.044162663 0.0140890950
## 1769469: -1.786385 0.009565163 0.0042702613
## 1769470: -1.786385 -0.003614836 0.0000000000
## 1769471: -1.786385 -0.003614836 0.0000000000
## 1769472: -1.786385 -0.003614836 0.0008876992# sst_lm <-
# sst[!is.na(Mean_Ampl) & !is.na(asst) & sd != 0,
# {
# a <- summary(lm(Mean_Ampl ~ asst))
# list(estimate = a$coefficients[2, 1],
# se = a$coefficients[2, 2])
# },
# by = .(lon, lat, month)]
# save(sst_lm, file = "sst_lm.rda")
load("sst_lm.rda")
# sst_ARMA <-
# sst[!is.na(Mean_Ampl) & !is.na(asst) & sd != 0,
# {
# a <- arima(Mean_Ampl, xreg = asst,
# order = c(1, 0, 1), method = "ML")
# list(estimate = coef(a)[[4]],
# se = sqrt(diag(a$var.coef))[4])
# },
# by = .(lon, lat, month)]
# save(sst_ARMA, file = "sst_ARMA.rda")
load("sst_ARMA.rda")
sst_ARMA[, model := "ARMA"]## lon lat month estimate se model
## 1: 165.00 -76.070 1 -0.5710952 6.160772 ARMA
## 2: 168.75 -76.070 1 -1.0979065 6.955729 ARMA
## 3: 172.50 -76.070 1 16.3959292 6.298035 ARMA
## 4: 176.25 -76.070 1 5.9687953 9.500008 ARMA
## 5: 180.00 -76.070 1 -1.4928169 7.130592 ARMA
## ---
## 35617: 333.75 87.159 12 -227.2418053 576.201620 ARMA
## 35618: 337.50 87.159 12 -209.5381702 184.037547 ARMA
## 35619: 341.25 87.159 12 -265.9381561 249.217011 ARMA
## 35620: 345.00 87.159 12 -464.5931503 675.701240 ARMA
## 35621: 356.25 87.159 12 -1662.0564847 2741.928059 ARMAsst_lm[, model := "lm"]## lon lat month estimate se model
## 1: 165.00 -76.070 1 3.465502 7.369890 lm
## 2: 168.75 -76.070 1 6.279632 7.834251 lm
## 3: 172.50 -76.070 1 9.603780 9.051243 lm
## 4: 176.25 -76.070 1 13.800397 9.930859 lm
## 5: 180.00 -76.070 1 7.110357 7.979328 lm
## ---
## 35617: 333.75 87.159 12 -418.695999 589.880315 lm
## 35618: 337.50 87.159 12 -66.482428 156.488025 lm
## 35619: 341.25 87.159 12 -67.993298 205.546154 lm
## 35620: 345.00 87.159 12 -288.471974 677.332404 lm
## 35621: 356.25 87.159 12 -1405.224116 3258.036442 lmsst_regr <- rbind(sst_ARMA, sst_lm)
remove(sst_ARMA, sst_lm)
sst_regr[, t := abs(estimate/se)]## lon lat month estimate se model t
## 1: 165.00 -76.070 1 -0.5710952 6.160772 ARMA 0.09269863
## 2: 168.75 -76.070 1 -1.0979065 6.955729 ARMA 0.15784205
## 3: 172.50 -76.070 1 16.3959292 6.298035 ARMA 2.60334024
## 4: 176.25 -76.070 1 5.9687953 9.500008 ARMA 0.62829372
## 5: 180.00 -76.070 1 -1.4928169 7.130592 ARMA 0.20935386
## ---
## 71238: 333.75 87.159 12 -418.6959993 589.880315 lm 0.70979822
## 71239: 337.50 87.159 12 -66.4824280 156.488025 lm 0.42484035
## 71240: 341.25 87.159 12 -67.9932976 205.546154 lm 0.33079333
## 71241: 345.00 87.159 12 -288.4719743 677.332404 lm 0.42589425
## 71242: 356.25 87.159 12 -1405.2241162 3258.036442 lm 0.43131013ggplot(sst_regr[model == "lm" & abs(estimate) < 40], aes(lon, lat)) +
geom_raster(aes(fill = estimate)) +
# geom_contour(aes(z = abs(estimate/se)), breaks = 2:100, color = "black") +
geom_point(size = 0.1, shape = 3, alpha = 0.4,
data = sst_regr[model == "lm" & t > 2]) +
world3 + coord_quickmap() +
theme_elio +
facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) +
scale_fill_gradient2(high = muted("red"), low = muted("blue"),
name = "Regresión (LM)", limits = c(-40, 40))
Si quiero ser un poco más riguroso y tener en cuenta la estructura de correlación que seguro existe en las series, puedo usar un modelo ARMA(1, 1) para estimar la pendiente.
ggplot(sst_regr[model == "ARMA" & abs(estimate) < 40], aes(lon, lat)) +
geom_raster(aes(fill = estimate)) +
# geom_contour(aes(z = abs(estimate/se)), breaks = 2:100, color = "black") +
geom_point(size = 0.1, shape = 3, alpha = 0.4,
data = sst_regr[model == "ARMA" & t > 2]) +
world3 + coord_quickmap() +
theme_elio +
facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) +
scale_fill_gradient2(high = muted("red"), low = muted("blue"),
name = "Regresión (ARMA)", limits = c(-40, 40))
En ambas figuras las cruces marcan regiones donde la magnitud de la regresión es mayor a 2 veces el error estandar de la misma (~ p-valor < 0.05).
(Alternativa: hacer test t-student para cada punto de latitud entre años activos e inactivos. )
# test <-
# sst[!is.na(Mean_Ampl) & !is.na(asst) & sd != 0,
# t.test(asst[Anom == "Activo"], y = asst[Anom == "Inactivo"]),
# by = .(lat, lon, month)]
#
# ggplot(test, aes(lon, lat)) +
# geom_raster(aes(fill = estimate)) +
# geom_point(aes(alpha = p.value < 0.05), size = 0.1, shape = 3) +
# scale_alpha_discrete(range = c(0, 0.7)) +
# scale_fill_gradient2(high = muted("red"), low = muted("blue")) +
# facet_wrap(~month, labeller = labeller(month = month.abb), ncol = 3) +
# coord_quickmap() +
# world3 +
# theme_elio